/* DISPLACEMENT AND STRAIN AT DEPTH (PART-B) DUE TO BURIED FINITE
   FAULT IN A SEMIINFINITE MEDIUM */

#include <cmath>
#include "../Iteration_Params.hxx"

Okada::Displacement Iteration_Params::ub
(const FTensor::Tensor1<double,3> &dislocation,
 const double &alp3,
 const double &sin_dip,
 const double &cos_dip,
 const double &cos_dip2,
 const double &sin_dip2,
 const double &sin_cos_dip) const
{
  const double rd(r+d), d11(1/(r*rd)), aj2(xi*y/rd*d11), aj5(-(d+y*y/rd)*d11);
  double ai3, ai4, ak1, ak3, aj3, aj6;
  if(cos_dip!=0)
    {
      if(xi==0)
        {
          ai4=0;
        }
      else
        {
          const double x(sqrt(xi2+q2));
          ai4=1/cos_dip2*( xi/rd*sin_cos_dip
                           +2*atan((et*(x+q*cos_dip)+x*(r+x)*sin_dip)/(xi*(r+x)*cos_dip)));
        }
      ai3=(y*cos_dip/rd-ale+sin_dip*log(rd))/cos_dip2;
      ak1=xi*(d11-y11*sin_dip)/cos_dip;
      ak3=(q*y11-y*d11)/cos_dip;
      aj3=(ak1-aj2*sin_dip)/cos_dip;
      aj6=(ak3-aj5*sin_dip)/cos_dip;
    }
  else
    {
      const double rd2(rd*rd);
      ai3=(et/rd+y*q/rd2-ale)/2;
      ai4=xi*y/rd2/2;
      ak1=xi*q/rd*d11;
      ak3=sin_dip/rd*(xi2*d11-1);
      aj3=-xi/rd2*(q2*d11-0.5);
      aj6=-y/rd2*(xi2*d11-0.5);
    }
  const double xy(xi*y11), ai1(-xi/rd*cos_dip-ai4*sin_dip),
    ai2(log(rd)+ai3*sin_dip), ak2(1/r+ak3*sin_dip),
    ak4(xy*cos_dip-ak1*sin_dip), aj1(aj5*cos_dip-aj6*sin_dip),
    aj4(-xy-aj2*cos_dip+aj3*sin_dip);

  Okada::Displacement result;
  FTensor::Index<'a',3> a;
  FTensor::Index<'b',3> b;
  result.first(a)=0;
  result.second(a,b)=0;

  const double qx(q*x11), qy(q*y11);

  FTensor::Tensor1<double,3> first;
  FTensor::Tensor2<double,3,3> second;
  const double pi(4*atan(1.0));
  /* STRIKE-SLIP CONTRIBUTION */
  if(dislocation(0)!=0)
    {
      first(0)=-xi*qy-tt -alp3*ai1*sin_dip;
      first(1)=-q/r      +alp3*y/rd*sin_dip;
      first(2)= q*qy     -alp3*ai2*sin_dip;
      second(0,0)= xi2*q*y32 -alp3*aj1*sin_dip;
      second(0,1)= xi*q/r3   -alp3*aj2*sin_dip;
      second(0,2)=-xi*q2*y32 -alp3*aj3*sin_dip;
      second(1,0)=-xi*fy-d*x11 +alp3*(xy+aj4)*sin_dip;
      second(1,1)=-ey          +alp3*(1/r+aj5)*sin_dip;
      second(1,2)= q*fy        -alp3*(qy-aj6)*sin_dip;
      second(2,0)=-xi*fz-y*x11 +alp3*ak1*sin_dip;
      second(2,1)=-ez          +alp3*y*d11*sin_dip;
      second(2,2)= q*fz        +alp3*ak2*sin_dip;

      result.first(a)=first(a)*dislocation(0)/(2*pi);
      result.second(a,b)=second(a,b)*dislocation(0)/(2*pi);
    }
  /* DIP-SLIP CONTRIBUTION */
  if(dislocation(1)!=0)
    {
      first(0)=-q/r      +alp3*ai3*sin_cos_dip;
      first(1)=-et*qx-tt -alp3*xi/rd*sin_cos_dip;
      first(2)= q*qx     +alp3*ai4*sin_cos_dip;
      second(0,0)= xi*q/r3     +alp3*aj4*sin_cos_dip;
      second(0,1)= et*q/r3+qy  +alp3*aj5*sin_cos_dip;
      second(0,2)=-q2/r3       +alp3*aj6*sin_cos_dip;
      second(1,0)=-ey          +alp3*aj1*sin_cos_dip;
      second(1,1)=-et*gy-xy*sin_dip +alp3*aj2*sin_cos_dip;
      second(1,2)= q*gy        +alp3*aj3*sin_cos_dip;
      second(2,0)=-ez          -alp3*ak3*sin_cos_dip;
      second(2,1)=-et*gz-xy*cos_dip -alp3*xi*d11*sin_cos_dip;
      second(2,2)= q*gz        -alp3*ak4*sin_cos_dip;

      result.first(a)+=first(a)*dislocation(1)/(2*pi);
      result.second(a,b)+=second(a,b)*dislocation(1)/(2*pi);
    }
  /* TENSILE-FAULT CONTRIBUTION */
  if(dislocation(2)!=0)
    {
      first(0)= q*qy           -alp3*ai3*sin_dip2;
      first(1)= q*qx           +alp3*xi/rd*sin_dip2;
      first(2)= et*qx+xi*qy-tt -alp3*ai4*sin_dip2;
      second(0,0)=-xi*q2*y32 -alp3*aj4*sin_dip2;
      second(0,1)=-q2/r3     -alp3*aj5*sin_dip2;
      second(0,2)= q*q2*y32  -alp3*aj6*sin_dip2;
      second(1,0)= q*fy -alp3*aj1*sin_dip2;
      second(1,1)= q*gy -alp3*aj2*sin_dip2;
      second(1,2)=-q*hy -alp3*aj3*sin_dip2;
      second(2,0)= q*fz +alp3*ak3*sin_dip2;
      second(2,1)= q*gz +alp3*xi*d11*sin_dip2;
      second(2,2)=-q*hz +alp3*ak4*sin_dip2;

      result.first(a)+=first(a)*dislocation(2)/(2*pi);
      result.second(a,b)+=second(a,b)*dislocation(2)/(2*pi);
    }
  return result;
}
